/*-
 * Copyright (c) 2014 Shang-Rong Cai. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * $ trigrid.c    Rev. 0.1.4    2014-05-30T11:03+08:00 $
 */

#include <stdlib.h>
#include <stdio.h>

#include <pregrid.h>
#include <trigrid.h>


int grid_triangle_create(int m, int n, int *vlst);

int grid_triangle_stdout(int lstz, int *vlst);

int grid_triangle_vtkout(FILE *fpof, int lstz, int *vlst);

int point_coordinate(double x0, double x1, int m, double y0, double y1, int n, double *coor);

int point_stdout(int lstz, double *coor);

int point_vtkout(FILE *fpof, int lstz, double *coor);


int count_vertex_cell_triangle(int vlstz, int *vlst, int vclz, int *vcz_lst);

int vcz_lst_stdout(size_t vclz, int *vcz_lst);

int vcz_lst_maximum(size_t vclz, int *vcz_lst);


int compare_index(const void *idxa, const void *idxb);

int count_vertex_degrees_neighbors(int vlstz, int *vlst, int vclz, int bkz, int *nbz_lst, int *nbr_lst);

int vertex_degrees_neighbors_stdout(int vclz, int bkz, int *nbz_lst, int *nbr_lst);


int out_degree_sum(int vclz, int *nbz_lst);

int cidk_structure(int vclz, int bkz, int *nbz_lst, int *nbr_lst, cidk_t *cidk);


int trigrid_create(double x0, double x1, double y0, double y1, size_t mz, size_t nz, trigrid_t *tgrd)
{
    int ierr = ERRNO_UNDEFINED;

    tgrd->pdim = TRIGRID_PROBLEM_DIMENSION;
    tgrd->sdim = TRIGRID_SPATIAL_DIMENSION;

    tgrd->rctr = TRIGRID_DEFAULT_RCTR;
    tgrd->trvz = TRIGRID_DEFAULT_TRVZ;

    tgrd->vtxz = (mz + 1) * (nz + 1);
    tgrd->elmz = mz * nz * tgrd->rctr;

    tgrd->coor = (double*) calloc(tgrd->vtxz * tgrd->sdim, sizeof(double));
    tgrd->vlst = (int*) calloc(tgrd->elmz * tgrd->trvz, sizeof(int));

    /* Generate Vertices Coordinates */
    ierr = point_coordinate(x0, x1, mz, y0, y1, nz, tgrd->coor);

    /* Generate Cell */
    ierr = grid_triangle_create(mz, nz, tgrd->vlst);

    return ierr;
}


int trigrid_delete(trigrid_t *tgrd)
{
    int ierr = ERRNO_UNDEFINED;

    tgrd->pdim = 0;
    tgrd->sdim = 0;

    tgrd->rctr = 0;
    tgrd->trvz = 0;

    free(tgrd->coor);
    free(tgrd->vlst);

    return ierr;
}


int trigrid_stdout(const trigrid_t *tgrd)
{
    int ierr = ERRNO_UNDEFINED;

    printf("Vertex List:\n");
    ierr = point_stdout(tgrd->vtxz, tgrd->coor);

    printf("Element List:\n");
    ierr = grid_triangle_stdout(tgrd->elmz, tgrd->vlst);

    return ierr;
}


int trigrid_save_vtk(const trigrid_t *tgrd, const char *vtk_file)
{
    int ierr = ERRNO_UNDEFINED;

    FILE *fpof;
    size_t zi;

    int cell_type;

    fpof = fopen(vtk_file, "w");

    fprintf(fpof, "# vtk DataFile Version 2.0\n");
    fprintf(fpof, "Simple Test Triangle Grid\n");
    fprintf(fpof, "ASCII\n");
    fprintf(fpof, "DATASET UNSTRUCTURED_GRID\n");

    fprintf(fpof, "POINTS %zd float\n", tgrd->vtxz);
    ierr = point_vtkout(fpof, tgrd->vtxz, tgrd->coor);

    fprintf(fpof, "CELLS %zd %zd\n", tgrd->elmz, tgrd->elmz * (tgrd->trvz + 1));
    ierr = grid_triangle_vtkout(fpof, tgrd->elmz, tgrd->vlst);

    fprintf(fpof, "CELL_TYPES %zd", tgrd->elmz);
    cell_type = 5;

    for (zi = 0; zi < tgrd->elmz; zi++)
    {
        if (zi % 20 == 0)
        {
            fprintf(fpof, "\n");
        }

        fprintf(fpof, "%4d", cell_type);
    }

    fclose(fpof);

    return ierr;
}


int trigrid_create_cidk(const trigrid_t *tgrd, cidk_t *cidk)
{
    int ierr = ERRNO_UNDEFINED;

    int *vcz_lst;  /* List number of cells for vertice belong to */
    int max_vcz;  /* Maximum number of cells of vertex belong to */

    int *nbz_lst;  /* Number of (closed) neighbors of vertices */
    int *nbr_lst;  /* (Closed) Neighbors of vertices */
    int deg_sum;  /* Degree sum */

    /* Count Vertex Belong To Cells */
    vcz_lst = (int*) calloc(tgrd->vtxz, sizeof(int));
    ierr = count_vertex_cell_triangle(tgrd->elmz, tgrd->vlst, tgrd->vtxz, vcz_lst);
    max_vcz = vcz_lst_maximum(tgrd->vtxz, vcz_lst);
    free(vcz_lst);

#ifdef DEBUG_INFO
    printf("Number of Cells To Each Vertices:\n");
    ierr = vcz_lst_stdout(tgrd->vtxz, vcz_lst);
    printf("Maximum: %d\n", max_vcz);
#endif

    /* Allocate Memory nbz_lst and nbr_lst */
    nbz_lst = (int*) calloc(tgrd->vtxz, sizeof(int));
    nbr_lst = (int*) calloc(tgrd->vtxz * tgrd->sdim * max_vcz, sizeof(int));

    /* Consider As Simple Di-Graph For Counting Degree */
    ierr = count_vertex_degrees_neighbors(tgrd->elmz, tgrd->vlst, tgrd->vtxz, tgrd->sdim * max_vcz, nbz_lst, nbr_lst);
#ifdef DEBUG_INFO
    ierr = vertex_degrees_neighbors_stdout(tgrd->vtxz, sdim * max_vcz, nbz_lst, nbr_lst);
#endif

    deg_sum = out_degree_sum(tgrd->vtxz, nbz_lst);
#ifdef DEBUG_INFO
    printf("Degree Sum: %d\n", deg_sum);
#endif

    /* Sparse Structure Data */
    ierr = cidk_init(tgrd->vtxz, tgrd->vtxz, (size_t) deg_sum, ROW_WISE, cidk);
    ierr = cidk_structure(tgrd->vtxz, tgrd->sdim * max_vcz, nbz_lst, nbr_lst, cidk);
#ifdef DEBUG_INFO
    printf("Sparse Matrix Value:\n");
    ierr = cidk_stdout(cidk);
#endif

    free(nbz_lst);
    free(nbr_lst);

    return ierr;
}


int trigrid_delete_cidk(cidk_t *cidk)
{
    return cidk_dest(cidk);;
}


int grid_triangle_create(int m, int n, int *vlst)
{
    int ierr = ERRNO_UNDEFINED;

    int ii, ij;
    int eidx, vidx;

    if (m < 1 || n < 1)
    {
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
        goto EXIT_GRID_TRIANGLE;
    }

    for (ii = 0; ii < m; ii++)
    {
        for (ij = 0; ij < n; ij++)
        {
            eidx = ii * n + ij;
            vidx = ii * (n + 1) + ij;

            vlst[eidx * 6] = vidx;
            vlst[eidx * 6 + 1] = vidx + 1;
            vlst[eidx * 6 + 2] = vidx + n + 1;

            vlst[eidx * 6 + 3] = vidx + 1;
            vlst[eidx * 6 + 4] = vidx + n + 2;
            vlst[eidx * 6 + 5] = vidx + n + 1;
        }
    }

    EXIT_GRID_TRIANGLE:

    return ierr;
}


int grid_triangle_stdout(int lstz, int *vlst)
{
    int ierr = ERRNO_UNDEFINED;

    int ii;

    for (ii = 0; ii < lstz; ii++)
    {
        printf("%7d: %7d %7d %7d\n", ii, vlst[ii * 3], vlst[ii * 3 + 1], vlst[ii * 3 + 2]);
    }

    return ierr;
}


int grid_triangle_vtkout(FILE *fpof, int lstz, int *vlst)
{
    int ierr = ERRNO_UNDEFINED;

    int ii;

    for (ii = 0; ii < lstz; ii++)
    {
        fprintf(fpof, "3 %7d %7d %7d\n", vlst[ii * 3], vlst[ii * 3 + 1], vlst[ii * 3 + 2]);
    }

    return ierr;
}


int point_coordinate(double x0, double x1, int m, double y0, double y1, int n, double *coor)
{
    int ierr = ERRNO_UNDEFINED;

    int ii, ij, idx;

    for (ii = 0; ii < m + 1; ii++)
    {
        for (ij = 0; ij < n + 1; ij++)
        {
            idx = ii * (n + 1) + ij;

            coor[idx * 3] = x0 + (x1 - x0) * (double) ii / (double) m;
            coor[idx * 3 + 1] = y0 + (y1 - y0) * (double) ij / (double) n;
            coor[idx * 3 + 2] = 0;
        }
    }

    return ierr;
}


int point_stdout(int lstz, double *coor)
{
    int ierr = ERRNO_UNDEFINED;

    int ii;

    for (ii = 0; ii < lstz; ii++)
    {
        printf("%7d: % .7lE % .7lE % .7lE\n", ii, coor[ii * 3], coor[ii * 3 + 1], coor[ii * 3 + 2]);
    }

    return ierr;
}


int point_vtkout(FILE *fpof, int lstz, double *coor)
{
    int ierr = ERRNO_UNDEFINED;

    int ii;

    for (ii = 0; ii < lstz; ii++)
    {
        fprintf(fpof, " % .7lE % .7lE % .7lE\n", coor[ii * 3], coor[ii * 3 + 1], coor[ii * 3 + 2]);
    }

    return ierr;
}


int count_vertex_cell_triangle(int vlstz, int *vlst, int vclz, int *vcz_lst)
{
    int ierr = ERRNO_UNDEFINED;

    int ii;

    for (ii = 0; ii < vclz; ii++)
    {
        vcz_lst[ii] = 0;
    }

    for (ii = 0; ii < vlstz; ii++)
    {
        vcz_lst[vlst[ii * 3]]++;
        vcz_lst[vlst[ii * 3 + 1]]++;
        vcz_lst[vlst[ii * 3 + 2]]++;
    }

    return ierr;
}


int vcz_lst_stdout(size_t vclz, int *vcz_lst)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;

    for (zi = 0; zi < vclz; zi++)
    {
        printf("%7zd: %3d\n", zi, vcz_lst[zi]);
    }

    return ierr;
}


int vcz_lst_maximum(size_t vclz, int *vcz_lst)
{
    int max_value;
    size_t zi;

    max_value = vcz_lst[0];

    for (zi = 1; zi < vclz; zi++)
    {
        if (max_value < vcz_lst[zi])
        {
            max_value = vcz_lst[zi];
        }
    }

    return max_value;
}


int compare_index(const void *idxa, const void *idxb)
{
    int index_a = *(const int*) idxa;
    int index_b = *(const int*) idxb;

    return ((index_a > index_b) - (index_a < index_b));
}


int count_vertex_degrees_neighbors(int vlstz, int *vlst, int vclz, int bkz, int *nbz_lst, int *nbr_lst)
{
    int ierr = ERRNO_UNDEFINED;

    int ii, ij;

    int idx_self, idx_nbr0, idx_nbr1;
    int buf_idx, cnt_cur;

    /* Initialize Zero For Degrees */
    for (ii = 0; ii < vclz; ii++)
    {
        nbz_lst[ii] = 0;
    }

    for (ii = 0; ii < vlstz; ii++)
    {
        /* Vertex 0 */
        idx_self = vlst[ii * 3];
        idx_nbr0 = vlst[ii * 3 + 1];
        idx_nbr1 = vlst[ii * 3 + 2];

        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_self;
        nbz_lst[idx_self]++;
        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_nbr0;
        nbz_lst[idx_self]++;
        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_nbr1;
        nbz_lst[idx_self]++;

        /* Vertex 1 */
        idx_self = vlst[ii * 3 + 1];
        idx_nbr0 = vlst[ii * 3];
        idx_nbr1 = vlst[ii * 3 + 2];

        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_self;
        nbz_lst[idx_self]++;
        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_nbr0;
        nbz_lst[idx_self]++;
        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_nbr1;
        nbz_lst[idx_self]++;

        /* Vertex 2 */
        idx_self = vlst[ii * 3 + 2];
        idx_nbr0 = vlst[ii * 3];
        idx_nbr1 = vlst[ii * 3 + 1];

        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_self;
        nbz_lst[idx_self]++;
        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_nbr0;
        nbz_lst[idx_self]++;
        nbr_lst[idx_self * bkz + nbz_lst[idx_self]] = idx_nbr1;
        nbz_lst[idx_self]++;
    }


    /* Sort And (Check Stdout) */
    for (ii = 0; ii < vclz; ii++)
    {
        qsort((void*) &nbr_lst[ii * bkz], (size_t) nbz_lst[ii], sizeof(int), compare_index);

    #ifdef DEBUG_INFO
        printf("%3d:", ii);
        for (ij = 0; ij < nbz_lst[ii]; ij++)
        {
            printf(" %3d", nbr_lst[ii * bkz + ij]);
        }
        printf("\n");
    #endif
    }

    /* Reduce */
    for (ii = 0; ii < vclz; ii++)
    {
        buf_idx = nbr_lst[ii * bkz];
        /*
        printf("%3d: %3d", ii, buf_idx);
        */
        cnt_cur = 1;

        for (ij = 1; ij < nbz_lst[ii]; ij++)
        {
            if (buf_idx != nbr_lst[ii * bkz + ij])
            {
                buf_idx = nbr_lst[ii * bkz + ij];
                nbr_lst[ii * bkz + cnt_cur] = buf_idx;
                cnt_cur++;

                /*
                printf(" %3d", buf_idx);
                */
            }
        }

        /*
        printf(" (%5d)\n", cnt_cur);
        */
        nbz_lst[ii] = cnt_cur;
    }

    return ierr;
}


int vertex_degrees_neighbors_stdout(int vclz, int bkz, int *nbz_lst, int *nbr_lst)
{
    int ierr = ERRNO_UNDEFINED;

    int ii, ij;

    for (ii = 0; ii < vclz; ii++)
    {
        printf("%3d: degree = %d;", ii, nbz_lst[ii]);

        for (ij = 0; ij < nbz_lst[ii]; ij++)
        {
            printf(" %3d", nbr_lst[ii * bkz + ij]);
        }

        printf("\n");
    }

    return ierr;
}


int out_degree_sum(int vclz, int *nbz_lst)
{
    int degsum = 0;

    int ii;

    for (ii = 0; ii < vclz; ii++)
    {
        degsum += nbz_lst[ii];
    }

    return degsum;
}


int cidk_structure(int vclz, int bkz, int *nbz_lst, int *nbr_lst, cidk_t *cidk)
{
    int ierr = ERRNO_UNDEFINED;

    int ii, ij;

    for (ii = 0; ii < vclz; ii++)
    {
        cidk->sptr[ii].head = cidk->e_sz;

        for (ij = 0; ij < nbz_lst[ii]; ij++)
        {
            cidk->item[cidk->e_sz].index = nbr_lst[ii * bkz + ij];
            cidk->item[cidk->e_sz].value = 0;
            cidk->e_sz++;
        }

        cidk->sptr[ii].tail = cidk->e_sz;
    }

    return ierr;
}
